#include <algorithm/analytical_model/CommonPara.h>
#include <algorithm/analytical_model/Constants.h>
#include <algorithm/analytical_model/HeliostatPara.h>
#include <algorithm/analytical_model/ScenePara.h>
#include <algorithm/analytical_model/SimulationPara.h>
#include <algorithm/analytical_model/algorithm/ModelPara.h>
#include <algorithm/analytical_model/algorithm/geometry/AllGeometry.h>
#include <algorithm/analytical_model/algorithm/geometry/Utils.h>
#include <algorithm/analytical_model/algorithm/simplified/HFLCALBasedCommon.h>
#include <algorithm/analytical_model/simplified/BHFLCALPara.h>
#include <algorithm/shared/data/scene/Heliostat.h>
#include <algorithm/shared/data/scene/Receiver.h>
#include <algorithm/shared/data/topology/geometry/Geometry.h>
#include <algorithm/shared/data/topology/geometry/Rectangle3D.h>
#include <utils/HasValue.h>
#include <utils/Math.h>
#include <utils/ParaUtils.h>
#include <utils/config.h>


#include <cmath>

using namespace solar;

template <>
void solar::computeSigma(SimulationPara* simulation_para, ScenePara* scene_para,
                         HeliostatPara* heliostat_para, CommonPara* common_para, Sun* sun,
                         Heliostat* heliostat, Receiver* receiver, BHFLCALPara* ptr)
{
    const auto [sigma_x, peak_flux_param] = internal::computeHFLCALBasedSigma(
        simulation_para, scene_para, heliostat_para, sun, heliostat, receiver);
    const auto geometry_heliostat = SOLAR_GET_VALUE(Heliostat, geometry_, *heliostat);
    const auto geometry_receiver = SOLAR_GET_VALUE(Receiver, geometry_, *receiver);
    const auto image_plane_normal =
        solar::getCenter(geometry_heliostat.get()) -
        solar::getFocusCenter(geometry_heliostat.get(), geometry_receiver.get());
    const auto heliostat_vertices = SOLAR_GET_VALUE(
        Rectangle3D, vertices_, *(dynamic_cast<Rectangle3D*>(geometry_heliostat.get())));
    const auto heliostat_length_vec = heliostat_vertices[3] - heliostat_vertices[0];
    const auto heliostat_width_vec = heliostat_vertices[1] - heliostat_vertices[0];
    const auto projected_heliostat_length =
        solar::projectTo(heliostat_length_vec, image_plane_normal).length();
    const auto projected_heliostat_width =
        solar::projectTo(heliostat_width_vec, image_plane_normal).length();

    const auto projected_length_ratio = projected_heliostat_length / projected_heliostat_width;
    auto sigma_y = sigma_x / (1 + std::log10(projected_length_ratio));
    if (sigma_y < 0)
    {
        sigma_y = sigma_x / const_expr::invSqrt(projected_length_ratio);
    }

    SOLAR_SET_VALUE(BHFLCALPara, sigma_x_, sigma_x, *ptr)
    SOLAR_SET_VALUE(BHFLCALPara, sigma_y_, sigma_y, *ptr)
    SOLAR_SET_VALUE(BHFLCALPara, flux_param_, peak_flux_param, *ptr)
}
